Klasyfikacja niezbalansowana, klasyfikatory zespołowe i wyjaśnialna AI

Wykorzystanie Google Colab

Jeśli korzystasz z Google Colab skopiuj plik feature_names.json do katalogu głównego projektu.

Ładowanie i eksploracja danych

Na tym laboratorium wykorzystamy zbiór danych Polish companies bankruptcy. Dotyczy on klasyfikacji, na podstawie danych z raportów finansowych, czy firma zbankrutuje w ciągu najbliższych kilku lat. Jest to zadanie szczególnie istotne dla banków, funduszy inwestycyjnych, firm ubezpieczeniowych itp., które z tego powodu zatrudniają licznie data scientistów. Zbiór zawiera 64 cechy, obliczone przez ekonomistów, którzy stworzyli ten zbiór, są one opisane na podlinkowanej wcześniej stronie. Dotyczą one zysków, posiadanych zasobów oraz długów firm.

Ściągnij i rozpakuj dane (Data Folder -> data.zip) do katalogu data obok tego notebooka. Znajduje się tam 5 plików w formacie .arff, wykorzystywanym głównie przez oprogramowanie Weka. Jest to program do "klikania" ML w interfejsie graficznym, jakiś czas temu popularny wśród mniej technicznych data scientistów. W Pythonie ładuje się je za pomocą bibliotek SciPy i Pandas.

Jeśli korzystasz z Linuksa możesz skorzystać z poniższych poleceń do pobrania i rozpakowania tych plików.

# !mkdir -p data
# !wget https://archive.ics.uci.edu/static/public/365/polish+companies+bankruptcy+data.zip -O data/data.zip

# !unzip data/data.zip -d data

W dalszej części laboratorium wykorzystamy plik 3year.arff, w którym na podstawie finansowych firmy po 3 latach monitorowania chcemy przewidywać, czy firma zbankrutuje w ciągu najbliższych 3 lat. Jest to dość realistyczny horyzont czasowy.

Dodatkowo w pliku feature_names.json znajdują się nazwy cech. Są bardzo długie, więc póki co nie będziemy z nich korzystać.

import json
import os

import numpy as np
import pandas as pd
from scipy.io import arff

data = arff.loadarff(os.path.join("data", "3year.arff"))

with open("feature_names.json") as file:
    feature_names = np.array(json.load(file))

X = pd.DataFrame(data[0])

Przyjrzyjmy się teraz naszym danym.

X.head()

	Attr1	Attr2	Attr3	Attr4	Attr5	Attr6	Attr7	Attr8	Attr9	Attr10	...	Attr56	Attr57	Attr58	Attr59	Attr60	Attr61	Attr62	Attr63	Attr64	class
0	0.174190	0.41299	0.14371	1.3480	-28.9820	0.60383	0.219460	1.1225	1.1961	0.46359	...	0.163960	0.375740	0.83604	0.000007	9.7145	6.2813	84.291	4.3303	4.0341	b'0'
1	0.146240	0.46038	0.28230	1.6294	2.5952	0.00000	0.171850	1.1721	1.6018	0.53962	...	0.027516	0.271000	0.90108	0.000000	5.9882	4.1103	102.190	3.5716	5.9500	b'0'
2	0.000595	0.22612	0.48839	3.1599	84.8740	0.19114	0.004572	2.9881	1.0077	0.67566	...	0.007639	0.000881	0.99236	0.000000	6.7742	3.7922	64.846	5.6287	4.4581	b'0'
3	0.024526	0.43236	0.27546	1.7833	-10.1050	0.56944	0.024526	1.3057	1.0509	0.56453	...	0.048398	0.043445	0.95160	0.142980	4.2286	5.0528	98.783	3.6950	3.4844	b'0'
4	0.188290	0.41504	0.34231	1.9279	-58.2740	0.00000	0.233580	1.4094	1.3393	0.58496	...	0.176480	0.321880	0.82635	0.073039	2.5912	7.0756	100.540	3.6303	4.6375	b'0'

5 rows × 65 columns

X.dtypes

Attr1     float64
Attr2     float64
Attr3     float64
Attr4     float64
Attr5     float64
           ...   
Attr61    float64
Attr62    float64
Attr63    float64
Attr64    float64
class      object
Length: 65, dtype: object

X.describe()

	Attr1	Attr2	Attr3	Attr4	Attr5	Attr6	Attr7	Attr8	Attr9	Attr10	...	Attr55	Attr56	Attr57	Attr58	Attr59	Attr60	Attr61	Attr62	Attr63	Attr64
count	10503.000000	10503.000000	10503.000000	10485.000000	1.047800e+04	10503.000000	10503.000000	10489.000000	10500.000000	10503.000000	...	1.050300e+04	10460.000000	10503.000000	10474.000000	10503.000000	9.911000e+03	10486.000000	1.046000e+04	10485.000000	10275.000000
mean	0.052844	0.619911	0.095490	9.980499	-1.347662e+03	-0.121159	0.065624	19.140113	1.819254	0.366093	...	6.638549e+03	-0.530082	-0.014817	3.848794	1.429319	5.713363e+02	13.935361	1.355370e+02	9.095149	35.766800
std	0.647797	6.427041	6.420056	523.691951	1.185806e+05	6.970625	0.651152	717.756745	7.581659	6.428603	...	5.989196e+04	55.978608	18.684047	190.201224	77.273270	3.715967e+04	83.704103	2.599116e+04	31.419096	428.298315
min	-17.692000	0.000000	-479.730000	0.002080	-1.190300e+07	-508.120000	-17.692000	-2.081800	-1.215700	-479.730000	...	-7.513800e+05	-5691.700000	-1667.300000	-198.690000	-172.070000	0.000000e+00	-6.590300	-2.336500e+06	-0.000156	-0.000102
25%	0.000686	0.253955	0.017461	1.040100	-5.207075e+01	0.000000	0.002118	0.431270	1.011275	0.297340	...	1.462100e+01	0.005137	0.006796	0.875560	0.000000	5.533150e+00	4.486075	4.073700e+01	3.062800	2.023350
50%	0.043034	0.464140	0.198560	1.605600	1.579300e+00	0.000000	0.050945	1.111000	1.199000	0.515500	...	8.822900e+02	0.051765	0.106880	0.953060	0.002976	9.952100e+00	6.677300	7.066400e+01	5.139200	4.059300
75%	0.123805	0.689330	0.419545	2.959500	5.608400e+01	0.072584	0.142275	2.857100	2.059100	0.725635	...	4.348900e+03	0.130010	0.271310	0.995927	0.240320	2.093600e+01	10.587500	1.182200e+02	8.882600	9.682750
max	52.652000	480.730000	17.708000	53433.000000	6.854400e+05	45.533000	52.652000	53432.000000	740.440000	11.837000	...	3.380500e+06	293.150000	552.640000	18118.000000	7617.300000	3.660200e+06	4470.400000	1.073500e+06	1974.500000	21499.000000

8 rows × 64 columns

feature_names

array(['net profit / total assets', 'total liabilities / total assets',
       'working capital / total assets',
       'current assets / short-term liabilities',
       '[(cash + short-term securities + receivables - short-term liabilities) / (operating expenses - depreciation)] * 365',
       'retained earnings / total assets', 'EBIT / total assets',
       'book value of equity / total liabilities', 'sales / total assets',
       'equity / total assets',
       '(gross profit + extraordinary items + financial expenses) / total assets',
       'gross profit / short-term liabilities',
       '(gross profit + depreciation) / sales',
       '(gross profit + interest) / total assets',
       '(total liabilities * 365) / (gross profit + depreciation)',
       '(gross profit + depreciation) / total liabilities',
       'total assets / total liabilities', 'gross profit / total assets',
       'gross profit / sales', '(inventory * 365) / sales',
       'sales (n) / sales (n-1)',
       'profit on operating activities / total assets',
       'net profit / sales', 'gross profit (in 3 years) / total assets',
       '(equity - share capital) / total assets',
       '(net profit + depreciation) / total liabilities',
       'profit on operating activities / financial expenses',
       'working capital / fixed assets', 'logarithm of total assets',
       '(total liabilities - cash) / sales',
       '(gross profit + interest) / sales',
       '(current liabilities * 365) / cost of products sold',
       'operating expenses / short-term liabilities',
       'operating expenses / total liabilities',
       'profit on sales / total assets', 'total sales / total assets',
       'constant capital / total assets', 'profit on sales / sales',
       '(current assets - inventory - receivables) / short-term liabilities',
       'total liabilities / ((profit on operating activities + depreciation) * (12/365))',
       'profit on operating activities / sales',
       'rotation receivables + inventory turnover in days',
       '(receivables * 365) / sales', 'net profit / inventory',
       '(current assets - inventory) / short-term liabilities',
       '(inventory * 365) / cost of products sold',
       'EBITDA (profit on operating activities - depreciation) / total assets',
       'EBITDA (profit on operating activities - depreciation) / sales',
       'current assets / total liabilities',
       'short-term liabilities / total assets',
       '(short-term liabilities * 365) / cost of products sold)',
       'equity / fixed assets', 'constant capital / fixed assets',
       'working capital', '(sales - cost of products sold) / sales',
       '(current assets - inventory - short-term liabilities) / (sales - gross profit - depreciation)',
       'total costs / total sales', 'long-term liabilities / equity',
       'sales / inventory', 'sales / receivables',
       '(short-term liabilities * 365) / sales',
       'sales / short-term liabilities', 'sales / fixed assets'],
      dtype='<U115')

DataFrame zawiera 64 atrybuty numeryczne o zróżnicowanych rozkładach wartości oraz kolumnę "class" typu bytes z klasami 0 i 1. Wiemy, że mamy do czynienia z klasyfikacją binarną - klasa 0 to brak bankructwa, klasa 1 to bankructwo w ciągu najbliższych 3 lat. Przyjrzyjmy się dokładniej naszym danym.

Zadanie 1 (0.5 punktu)

1. Wyodrębnij klasy jako osobną zmienną typu pd.Series, usuwając je z macierzy X. Przekonwertuj go na liczby całkowite.
2. Narysuj wykres słupkowy (bar plot) częstotliwości obu klas w całym zbiorze. Upewnij się, że na osi X są numery lub nazwy klas, a oś Y ma wartości w procentach.

Uwaga: sugerowane jest użycie if w podpunkcie 1, żeby można było tę komórkę bezpiecznie odpalić kilka razy.

import matplotlib.pyplot as plt


# function for plotting class distribution
def plot_class_distribution(y):
    counts = y.value_counts()
    percentages = y.value_counts(normalize=True) * 100

    ax = percentages.plot.bar()

    ax.set_title("Classes Frequency")
    ax.set_xlabel("Class")
    ax.set_ylabel("Frequency [%]")

    labels = [f"{count} ({perc:.2f}%)" for count, perc in zip(counts, percentages)]
    ax.bar_label(ax.containers[0], labels=labels)
    plt.show()

# your_code

if "class" in X.columns:
    y = X.pop("class")

y = y.astype(int)

plot_class_distribution(y)

assert "class" not in X.columns

print("Solution is correct!")

Solution is correct!

Jak widać, klasa pozytywna jest w znacznej mniejszości, stanowi poniżej 5% zbioru. Taki problem nazywamy klasyfikacją niezbalansowaną (imbalanced classification). Mamy tu klasę dominującą (majority class) oraz klasę mniejszościową (minority class). Pechowo prawie zawsze interesuje nas ta druga, bo klasa większościowa jest trywialna. Przykładowo, 99% badanych jest zdrowych, a 1% ma niewykryty nowotwór - z oczywistych przyczyn chcemy wykrywać właśnie sytuację rzadką (problem diagnozy jako klasyfikacji jest zasadniczo zawsze niezbalansowany). W dalszej części laboratorium poznamy szereg konsekwencji tego zjawiska i metody na radzenie sobie z nim.

Mamy sporo cech, wszystkie numeryczne. Ciekawe, czy mają wartości brakujące, a jeśli tak, to ile. Policzymy to z pomocą biblioteki Pandas i metody .isna(). Domyślnie operuje ona na kolumnach, jak większość metod w Pandasie. Sumę wartości per kolumna zwróci nam metoda .sum(). Jeżeli podzielimy to przez liczbę wierszy len(X), to otrzymamy ułamek wartości brakujących w każdej kolumnie.

Pandas potrafi też stworzyć wykres, z pomocą funkcji np. .plot.hist() czy .plot.bar(). Przyjmują one opcje formatowania wykresu, z których korzysta pod spodem biblioteka matplotlib.

na_perc = X.isna().sum() / len(X)
na_perc.plot.bar(title="Fraction of missing values per column", figsize=(15, 5))

<Axes: title={'center': 'Fraction of missing values per column'}>

Jak widać, cecha 37 ma bardzo dużo wartości brakujących, podczas gdy pozostałe cechy mają raczej niewielką ich liczbę. W takiej sytuacji najlepiej usunąć tę cechę, a pozostałe wartości brakujące uzupełnić / imputować (impute). Typowo wykorzystuje się do tego wartość średnią lub medianę z danej kolumny. Ale uwaga - imputacji dokonuje się dopiero po podziale na zbiór treningowy i testowy! W przeciwnym wypadku wykorzystywalibyśmy dane ze zbioru testowego, co sztucznie zawyżyłoby wyniki. Jest to błąd metodologiczny - wyciek danych (data leakage).

Podział na zbiór treningowy i testowy to pierwszy moment, kiedy niezbalansowanie danych nam przeszkadza. Jeżeli zrobimy to czysto losowo, to są spore szanse, że w zbiorze testowym będzie tylko klasa negatywna - w końcu jest jej aż >95%. Dlatego wykorzystuje się próbkowanie ze stratyfikacją (stratified sampling), dzięki któremu proporcje klas w zbiorze przed podziałem oraz obu zbiorach po podziale są takie same.

Zadanie 2 (0.75 punktu)

1. Usuń kolumnę "Attr37" ze zbioru danych.
2. Dokonaj podziału zbioru na treningowy i testowy w proporcjach 80%-20%, z przemieszaniem (shuffle), ze stratyfikacją, wykorzystując funkcję train_test_split ze Scikit-learn'a.
3. Uzupełnij wartości brakujące średnią wartością cechy z pomocą klasy SimpleImputer.

Uwaga:

• jak wcześniej, sugerowane jest użycie if w podpunkcie 1
• pamiętaj o uwzględnieniu stałego random_state=0, aby wyniki były reprodukowalne (reproducible)
• stratify oczekuje wektora klas
• wartości do imputacji trzeba wyestymować na zbiorze treningowym (.fit()), a potem zastosować te nauczone wartości na obu podzbiorach (treningowym i testowym)

# your_code

from sklearn.impute import SimpleImputer
from sklearn.model_selection import train_test_split

if "Attr37" in X:
    X = X.drop("Attr37", axis="columns")

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=0, shuffle=True, stratify=y
)

mean_imputer = SimpleImputer(strategy="mean")

X_train = mean_imputer.fit_transform(X_train)
X_test = mean_imputer.transform(X_test)

import numpy as np

assert "Attr37" not in X.columns
assert not np.any(np.isnan(X_train))
assert not np.any(np.isnan(X_test))

print("Solution is correct!")

Solution is correct!

Prosta klasyfikacja

Zanim przejdzie się do modeli bardziej złożonych, trzeba najpierw wypróbować coś prostego, żeby mieć punkt odniesienia. Tworzy się dlatego modele bazowe (baselines).

W naszym przypadku będzie to drzewo decyzyjne (decision tree). Jest to drzewo binarne z decyzjami if-else, prowadzącymi do klasyfikacji danego przykładu w liściu. Każdy podział w drzewie to pytanie postaci "Czy wartość cechy X jest większa lub równa Y?". Trening takiego drzewa to prosty algorytm zachłanny, bardzo przypomina budowę zwykłego drzewa binarnego. W każdym węźle wykonujemy:

1. Sprawdź po kolei wszystkie możliwe punkty podziału, czyli każdą (unikalną) wartość każdej cechy, po kolei.
2. Dla każdego przypadku podziel zbiór na 2 kawałki: niespełniający warunku (lewe dziecko) i spełniający warunek (prawe dziecko).
3. Oblicz jakość podziału według pewnej wybranej funkcji jakości. Im lepiej nasz if/else rozdziela klasy od siebie (im "czystsze" są węzły-dzieci), tym wyższa jakość. Innymi słowy, chcemy, żeby do jednego dziecka poszła jedna klasa, a do drugiego druga.
4. Wybierz podział o najwyższej jakości.

Taki algorytm wykonuje się rekurencyjnie, aż otrzymamy węzeł czysty (pure leaf), czyli taki, w którym są przykłady z tylko jednej klasy. Typowo wykorzystywaną funkcją jakości (kryterium podziału) jest entropia Shannona - im niższa entropia, tym bardziej jednolite są klasy w węźle (czyli wybieramy podział o najniższej entropii).

Powyższe wytłumaczenie algorytmu jest oczywiście nieformalne i dość skrótowe. Doskonałe tłumaczenie, z interaktywnymi wizualizacjami, dostępne jest tutaj. W formie filmów - tutaj oraz tutaj. Dla drzew do regresji - ten film.

Warto zauważyć, że taka konstrukcja prowadzi zawsze do overfittingu. Otrzymanie liści czystych oznacza, że mamy 100% dokładności na zbiorze treningowym, czyli perfekcyjnie przeuczony klasyfikator. W związku z tym nasze predykcje mają bardzo niski bias, ale bardzo dużą wariancję. Pomimo tego drzewa potrafią dać bardzo przyzwoite wyniki, a w celu ich poprawy można je regularyzować, aby mieć mniej "rozrośnięte" drzewo. Film dla zainteresowanych.

W tym wypadku AI to naprawdę tylko zbiór if'ów ;)

Mając wytrenowany klasyfikator, trzeba oczywiście sprawdzić, jak dobrze on sobie radzi. Tu natrafiamy na kolejny problem z klasyfikacją niezbalansowaną - zwykła celność (accuracy) na pewno nie zadziała! Typowo wykorzystuje się AUC, nazywane też AUROC (Area Under Receiver Operating Characteristic), bo metryka ta "widzi" i uwzględnia niezbalansowanie klas. Wymaga ona przekazania prawdopodobieństwa klasy pozytywnej, a nie tylko binarnej decyzji.

Bardzo dobre i bardziej szczegółowe wytłumaczenie, z interktywnymi wizualizacjami, można znaleć tutaj. Dla preferujących filmy - tutaj.

Co ważne, z definicji AUROC, trzeba tam użyć prawdopodobieństw klasy pozytywnej (klasy 1). W Scikit-learn'ie zwraca je metoda .predict_proba(), która w kolejnych kolumnach zwraca prawdopodobieństwa poszczególnych klas.

Zadanie 3 (0.75 punktu)

1. Wytrenuj klasyfikator drzewa decyzyjnego (klasa DecisionTreeClassifier). Użyj entropii jako kryterium podziału.
2. Oblicz i wypisz AUROC na zbiorze testowym dla drzewa decyzyjnego (funkcja roc_auc_score).
3. Skomentuj wynik - czy twoim zdaniem osiągnięty AUROC to dużo czy mało, biorąc pod uwagę możliwy zakres wartości tej metryki?

Uwaga:

• pamiętaj o użyciu stałego random_state=0
• jeżeli drzewo nie wyświetli się samo, użyj plt.show() z Matplotliba
• pamiętaj o tym, żeby przekazać do metryki AUROC prawdopodobieństwa klasy pozytywnej, a nie binarne predykcje!

from sklearn.metrics import roc_auc_score


# function for AUROC calculation
def get_auroc_score(model, X, y):

    y_score = model.predict_proba(X)[:, 1]

    auroc_score = roc_auc_score(y, y_score)

    print(f"AUROC score: {auroc_score}")

    return auroc_score

# your_code
from sklearn import tree
from sklearn.metrics import roc_auc_score
from sklearn.tree import DecisionTreeClassifier

model = DecisionTreeClassifier(random_state=0, criterion="entropy")

model.fit(X_train, y_train)

auroc = get_auroc_score(model, X_test, y_test)

AUROC score: 0.7266899766899767

fig = plt.figure(figsize=(200, 160))
_ = tree.plot_tree(model, filled=True, feature_names=X.columns, rounded=True)

assert auroc > 0.7

print("Solution is correct!")

Solution is correct!

// skomentuj tutaj

Moim zdaniem osiągnięty AUROC nie jest wysoki. Wytrenoway model ma "jakąś" moc klasyfikacyjną - AUROC > 0.5, więc nie jest to totalna losowość. Przejrzałem parę artykułów (1, 2, 3) i wyszło, że nie ma tak naprawdę czegoś takiego jak obiektywnie dobry AUROC, bo jest to zależne od sytuacji, ale często przyjmuje się, że ~0.7 to akceptowalny wynik - czyli na pewno da się lepiej.

Uczenie zespołowe, bagging, lasy losowe

Bardzo często wiele klasyfikatorów działających razem daje lepsze wyniki niż pojedynczy klasyfikator. Takie podejście nazywa się uczeniem zespołowym (ensemble learning). Istnieje wiele różnych podejść do tworzenia takich klasyfikatorów złożonych (ensemble classifiers).

Podstawową metodą jest bagging:

1. Wylosuj N (np. 100, 500, ...) próbek boostrapowych (boostrap sample) ze zbioru treningowego. Próbka boostrapowa to po prostu losowanie ze zwracaniem, gdzie dla wejściowego zbioru z M wierszami losujemy M próbek. Będą tam powtórzenia, średnio nawet 1/3, ale się tym nie przejmujemy.
2. Wytrenuj klasyfikator bazowy (base classifier) na każdej z próbek boostrapowych.
3. Stwórz klasyfikator złożony poprzez uśrednienie predykcji każdego z klasyfikatorów bazowych.

Typowo klasyfikatory bazowe są bardzo proste, żeby można było szybko wytrenować ich dużą liczbę. Prawie zawsze używa się do tego drzew decyzyjnych. Dla klasyfikacji uśrednienie wyników polega na głosowaniu - dla nowej próbki każdy klasyfikator bazowy ją klasyfikuje, sumuje się głosy na każdą klasę i zwraca najbardziej popularną decyzję.

Taki sposób ensemblingu zmniejsza wariancję klasyfikatora. Intuicyjnie, skoro coś uśredniamy, to siłą rzeczy będzie mniej rozrzucone, bo dużo ciężej będzie osiągnąć jakąś skrajność. Redukuje to też overfitting.

Lasy losowe (Random Forests) to ulepszenie baggingu. Zaobserwowano, że pomimo losowania próbek boostrapowych, w baggingu poszczególne drzewa są do siebie bardzo podobne (są skorelowane), używają podobnych cech ze zbioru. My natomiast chcemy zróżnicowania, żeby mieć niski bias - redukcją wariancji zajmuje się uśrednianie. Dlatego używa się metody losowej podprzestrzeni (random subspace method) - przy każdym podziale drzewa losuje się tylko pewien podzbiór cech, których możemy użyć do tego podziału. Typowo jest to pierwiastek kwadratowy z ogólnej liczby cech.

Zarówno bagging, jak i lasy losowe mają dodatkowo bardzo przyjemną własność - są mało czułe na hiperparametry, szczególnie na liczbę drzew. W praktyce wystarczy ustawić 500 czy 1000 drzew i będzie dobrze działać. Dalsze dostrajanie hiperparametrów może jeszcze trochę poprawić wyniki, ale nie tak bardzo, jak przy innych klasyfikatorach. Jest to zatem doskonały wybór domyślny, kiedy nie wiemy, jakiego klasyfikatora użyć.

Dodatkowo jest to problem embarassingly parallel - drzewa można trenować w 100% równolegle, dzięki czemu jest to dodatkowo wydajna obliczeniowo metoda.

Głębsze wytłumaczenie, z interaktywnymi wizualizacjami, można znaleźć tutaj. Dobrze tłumaczy je też ta seria filmów.

Zadanie 4 (0.5 punktu)

1. Wytrenuj klasyfikator Random Forest (klasa RandomForestClassifier). Użyj 500 drzew i entropii jako kryterium podziału.
2. Sprawdź AUROC na zbiorze testowym.
3. Skomentuj wynik w odniesieniu do drzewa decyzyjnego.

Uwaga: pamiętaj o ustawieniu random_state=0. Dla przyspieszenia ustaw n_jobs=-1 (użyje tylu procesów, ile masz dostępnych rdzeni procesora). Pamiętaj też o przekazaniu prawdopodobieństw do metryki AUROC.

# your_code
from sklearn.ensemble import RandomForestClassifier

model = RandomForestClassifier(
    n_estimators=500, criterion="entropy", random_state=0, n_jobs=-1
)

model.fit(X_train, y_train)

auroc = get_auroc_score(model, X_test, y_test)

AUROC score: 0.8994111948657404

assert auroc > 0.85

print("Solution is correct!")

Solution is correct!

// skomentuj tutaj

Otrzymaliśmy zdecydowanie lepszy wynik niż przy drzewie decyzyjnym. Bazując na ogólnie przyjętych normach jest to dobry wynik.

Jak zobaczymy poniżej, wynik ten możemy jednak jeszcze ulepszyć!

Oversampling, SMOTE

W przypadku zbiorów niezbalansowanych można dokonać balansowania (balancing) zbioru. Są tutaj 2 metody:

• undersampling: usunięcie przykładów z klasy dominującej
• oversampling: wygenerowanie dodatkowych przykładów z klasy mniejszościowej

Undersampling działa dobrze, kiedy niezbalansowanie jest niewielkie, a zbiór jest duży (możemy sobie pozwolić na usunięcie jego części). Oversampling typowo daje lepsze wyniki, istnieją dla niego bardzo efektywne algorytmy. W przypadku bardzo dużego niezbalansowania można zrobić oba.

Typowym algorytmem oversamplingu jest SMOTE (Synthetic Minority Oversampling TEchnique). Działa on następująco:

1. Idź po kolei po przykładach z klasy mniejszościowej
2. Znajdź k najbliższych przykładów dla próbki, typowo k=5
3. Wylosuj tylu sąsiadów, ile trzeba do oversamplingu, np. jeżeli chcemy zwiększyć klasę mniejszościową 3 razy (o 200%), to wylosuj 2 z 5 sąsiadów
4. Dla każdego z wylosowanych sąsiadów wylosuj punkt na linii prostej między próbką a tym sąsiadem. Dodaj ten punkt jako nową próbkę do zbioru

Taka technika generuje przykłady bardzo podobne do prawdziwych, więc nie zaburza zbioru, a jednocześnie pomaga klasyfikatorom, bo "zagęszcza" przestrzeń, w której znajduje się klasa pozytywna.

Algorytm SMOTE, jego warianty i inne algorytmy dla problemów niezbalansowanych implementuje biblioteka Imbalanced-learn.

Zadanie 5 (1 punkt)

Użyj SMOTE do zbalansowania zbioru treningowego (nie używa się go na zbiorze testowym!) (klasa SMOTE). Wytrenuj drzewo decyzyjne oraz las losowy na zbalansowanym zbiorze, użyj tych samych argumentów co wcześniej. Pamiętaj o użyciu wszędzie stałego random_state=0 oraz przekazaniu prawdopodobieństw do AUROC. Skomentuj wynik.

Wartość ROC drzewa decyzyjnego przypisz do zmiennej tree_roc, a lasu do forest_roc.

# your_code

from imblearn.over_sampling import SMOTE

smote = SMOTE(random_state=0)

X_reb, y_reb = smote.fit_resample(X_train, y_train)

plot_class_distribution(y_reb)

print(f"Decison tree:")
decision_tree = DecisionTreeClassifier(random_state=0, criterion="entropy")

decision_tree.fit(X_reb, y_reb)

tree_roc = get_auroc_score(decision_tree, X_test, y_test)

print()

print(f"Random forest:")

random_forest = RandomForestClassifier(
    n_estimators=500, criterion="entropy", random_state=0, n_jobs=-1
)

random_forest.fit(X_reb, y_reb)

forest_roc = get_auroc_score(random_forest, X_test, y_test)

Decison tree:
AUROC score: 0.70995670995671

Random forest:
AUROC score: 0.9047644274917003

assert 0.6 < tree_roc < 0.8
assert 0.8 < forest_roc < 0.95

// skomentuj tutaj

Wyniki są porównywalne do tych, które otrzymaliśmy poprzednio. AUROC można interpretować w ten sposób, że jest to prawdopodobieństwo, na to że spośród losowo wybranego TP i losowo wybranego TN, model zaklasyfikuje wyżej TP niż TN. Ta metryka nie hierarchizuje klas, przez co dobrze radzi sobie w przypadku ich niezbalansowania. Z tego powodu na zbiorze zrównoważonym, gdzie też nie uznaje jednej klasy za ważniejszej, radzi sobie podobnie.

W dalszej części laboratorium używaj zbioru po zastosowaniu SMOTE do treningu klasyfikatorów.

Dostrajanie (tuning) hiperparametrów

Lasy losowe są stosunkowo mało czułe na dobór hiperparametrów - i dobrze, bo mają ich dość dużo. Można zawsze jednak spróbować to zrobić, a w szczególności najważniejszy jest parametr max_features, oznaczający, ile cech losować przy każdym podziale drzewa. Typowo sprawdza się wartości z zakresu [0.1, 0.5].

W kwestii szybkości, kiedy dostrajamy hiperparametry, to mniej oczywiste jest, jakiego n_jobs użyć. Z jednej strony klasyfikator może być trenowany na wielu procesach, a z drugiej można trenować wiele klasyfikatorów na różnych zestawach hiperparametrów równolegle. Jeżeli nasz klasyfikator bardzo dobrze się uwspółbieżnia (jak Random Forest), to można dać mu nawet wszystkie rdzenie, a za to wypróbowywać kolejne zestawy hiperparametrów sekwencyjnie. Warto ustawić parametr verbose na 2 lub więcej, żeby dostać logi podczas długiego treningu i mierzyć czas wykonania. W praktyce ustawia się to metodą prób i błędów.

Zadanie 6 (1 punkt)

1. Dobierz wartość hiperparametru max_features:
   • użyj grid search z 5 foldami
   • wypróbuj wartości [0.1, 0.2, 0.3, 0.4, 0.5]
   • wybierz model o najwyższym AUROC (argument scoring)
2. Sprawdź, jaka była optymalna wartość max_features. Jest to atrybut wytrenowanego GridSearchCV.
3. Skomentuj wynik. Czy warto było poświęcić czas i zasoby na tę procedurę?
4. Wynik przypisz do zmiennej auroc.

Uwaga:

• pamiętaj, żeby jako estymatora przekazanego do grid search'a użyć instancji Random Forest, która ma już ustawione random_state=0 i n_jobs

# your_code

from sklearn.model_selection import GridSearchCV

parameters = {"max_features": [0.1, 0.2, 0.3, 0.4, 0.5]}

random_forest_estimator = RandomForestClassifier(
    n_estimators=500, criterion="entropy", random_state=0, n_jobs=-1
)

random_forest_GS = GridSearchCV(
    estimator=random_forest_estimator,
    cv=5,
    param_grid=parameters,
    verbose=2,
    scoring="roc_auc",
    n_jobs=-1,
)

random_forest_GS.fit(X_train, y_train)

Fitting 5 folds for each of 5 candidates, totalling 25 fits

GridSearchCV(cv=5,
             estimator=RandomForestClassifier(criterion='entropy',
                                              n_estimators=500, n_jobs=-1,
                                              random_state=0),
             n_jobs=-1, param_grid={'max_features': [0.1, 0.2, 0.3, 0.4, 0.5]},
             scoring='roc_auc', verbose=2)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

print(f'Optimal max_features value: {random_forest_GS.best_params_["max_features"]}')

auroc = get_auroc_score(random_forest_GS, X_test, y_test)

Optimal max_features value: 0.4
AUROC score: 0.9217020353383989

assert 0.9 <= auroc <= 0.95

print("Solution is correct!")

Solution is correct!

// skomentuj tutaj

Samo dostrajanie w moim przypadku nie zajęło zbyt dużo czasu - testowałem na dwóch maszynach i wyszło między 3-7 minut, więc akceptowalnie. Wynik się poprawił, więc jak najbardziej było warto.

W praktycznych zastosowaniach data scientist wedle własnego uznana, doświadczenia, dostępnego czasu i zasobów wybiera, czy dostrajać hiperparametry i w jak szerokim zakresie. Dla Random Forest na szczęście często może nie być znaczącej potrzeby, i za to go lubimy :)

Random Forest - podsumowanie

1. Model oparty o uczenie zespołowe
2. Kluczowe elementy:
   • bagging: uczenie wielu klasyfikatorów na próbkach boostrapowych
   • metoda losowej podprzestrzeni: losujemy podzbiór cech do każdego podziału drzewa
   • uśredniamy głosy klasyfikatorów
3. Dość odporny na overfitting, zmniejsza wariancję błędu dzięki uśrednianiu
4. Mało czuły na hiperparametry
5. Przeciętnie bardzo dobre wyniki, doskonały wybór domyślny przy wybieraniu algorytmu klasyfikacji

Boosting

Drugą bardzo ważną grupą algorytmów ensemblingu jest boosting, też oparty o drzewa decyzyjne. O ile Random Forest trenował wszystkie klasyfikatory bazowe równolegle i je uśredniał, o tyle boosting robi to sekwencyjnie. Drzewa te uczą się na całym zbiorze, nie na próbkach boostrapowych. Idea jest następująca: trenujemy drzewo decyzyjne, radzi sobie przeciętnie i popełnia błędy na częsci przykładów treningowych. Dokładamy kolejne, ale znające błędy swojego poprzednika, dzięki czemu może to uwzględnić i je poprawić. W związku z tym "boostuje" się dzięki wiedzy od poprzednika. Dokładamy kolejne drzewa zgodnie z tą samą zasadą.

Jak uczyć się na błędach poprzednika? Jest to pewna funkcja kosztu (błędu), którą chcemy zminimalizować. Zakłada się jakąś jej konkretną postać, np. squared error dla regresji, albo logistic loss dla klasyfikacji. Później wykorzystuje się spadek wzdłuż gradientu (gradient descent), aby nauczyć się, w jakim kierunku powinny optymalizować kolejne drzewa, żeby zminimalizować błędy poprzednika. Jest to konkretnie gradient boosting, absolutnie najpopularniejsza forma boostingu, i jeden z najpopularniejszych i osiągających najlepsze wyniki algorytmów ML.

Tyle co do intuicji. Ogólny algorytm gradient boostingu jest trochę bardziej skomplikowany. Bardzo dobrze i krok po kroku tłumaczy go ta seria filmów na YT. Szczególnie ważne implementacje gradient boostingu to XGBoost (Extreme Gradient Boosting) oraz LightGBM (Light Gradient Boosting Machine). XGBoost był prawdziwym przełomem w ML, uzyskując doskonałe wyniki i bardzo dobrze się skalując - był wykorzystany w CERNie do wykrywania cząstki Higgsa w zbiorze z pomiarów LHC mającym 10 milionów próbek. Jego implementacja jest dość złożona, ale dobrze tłumaczy ją inna seria filmików na YT.

Obecnie najczęściej wykorzystuje się LightGBM. Został stworzony przez Microsoft na podstawie doświadczeń z XGBoostem. Został jeszcze bardziej ulepszony i przyspieszony, ale różnice są głównie implementacyjne. Różnice dobrze tłumaczy ta prezentacja z konferencji PyData oraz prezentacja Microsoftu. Dla zainteresowanych - praktyczne aspekty LightGBM.

Zadanie 7 (0.5 punktu)

1. Wytrenuj klasyfikator LightGBM (klasa LGBMClassifier). Przekaż importance_type="gain" - przyda nam się to za chwilę.
2. Sprawdź AUROC na zbiorze testowym.
3. Skomentuj wynik w odniesieniu do wcześniejszych algorytmów.

Pamiętaj o random_state, n_jobs i prawdopodobieństwach dla AUROC.

# your_code

from lightgbm import LGBMClassifier

LGBM_clf = LGBMClassifier(importance_type="gain", random_state=0, n_jobs=-1)

LGBM_clf.fit(X_reb, y_reb)

auroc = get_auroc_score(LGBM_clf, X_test, y_test)

[LightGBM] [Info] Number of positive: 8006, number of negative: 8006
[LightGBM] [Info] Auto-choosing col-wise multi-threading, the overhead of testing was 0.003486 seconds.
You can set `force_col_wise=true` to remove the overhead.
[LightGBM] [Info] Total Bins 16065
[LightGBM] [Info] Number of data points in the train set: 16012, number of used features: 63
[LightGBM] [Info] [binary:BoostFromScore]: pavg=0.500000 -> initscore=0.000000
AUROC score: 0.9433748070111706

assert 0.9 <= auroc <= 0.97

print("Solution is correct!")

Solution is correct!

// skomentuj tutaj

LGBM osiągnał znakomity wynik, a ponadto, porównując z innymi klasyfikatorami, okazał się niesamowicie szybki, co tym bardziej działa na plus.

Boosting dzięki uczeniu na poprzednich drzewach redukuje nie tylko wariancję, ale też bias w błędzie, dzięki czemu może w wielu przypadkach osiągnąć lepsze rezultaty od lasu losowego. Do tego dzięki znakomitej implementacji LightGBM jest szybszy.

Boosting jest jednak o wiele bardziej czuły na hiperparametry niż Random Forest. W szczególności bardzo łatwo go przeuczyć, a większość hiperparametrów, których jest dużo, wiąże się z regularyzacją modelu. To, że teraz poszło nam lepiej z domyślnymi, jest rzadkim przypadkiem.

W związku z tym, że przestrzeń hiperparametrów jest duża, przeszukanie wszystkich kombinacji nie wchodzi w grę. Zamiast tego można wylosować zadaną liczbę zestawów hiperparametrów i tylko je sprawdzić - chociaż im więcej, tym lepsze wyniki powinniśmy dostać. Służy do tego RandomizedSearchCV. Co więcej, klasa ta potrafi próbkować rozkłady prawdopodobieństwa, a nie tylko sztywne listy wartości, co jest bardzo przydatne przy parametrach ciągłych.

Hiperparametry LightGBMa są dobrze opisane w oficjalnej dokumentacji: wersja krótsza i wersja dłuższa. Jest ich dużo, więc nie będziemy ich tutaj omawiać. Jeżeli chodzi o ich dostrajanie w praktyce, to przydatny jest oficjalny guide oraz dyskusje na Kaggle.

Zadanie 8 (1.5 punktu)

1. Zaimplementuj random search dla LightGBMa (klasa RandomizedSearchCV):
   • użyj tylu prób, na ile pozwalają twoje zasoby obliczeniowe, ale przynajmniej 30
   • przeszukaj przestrzeń hiperparametrów:

     param_grid = {
        "n_estimators": [100, 250, 500],
        "learning_rate": [0.05, 0.1, 0.2],
        "num_leaves": [31, 48, 64],
        "colsample_bytree": [0.8, 0.9, 1.0],
        "subsample": [0.8, 0.9, 1.0],
     }

2. Wypisz znalezione optymalne hiperparametry.
3. Wypisz raporty z klasyfikacji (funkcja classification_report), dla modelu LightGBM bez i z dostrajaniem hiperparametrów.
4. Skomentuj różnicę precyzji (precision) i czułości (recall) między modelami bez i z dostrajaniem hiperparametrów. Czy jest to pożądane zjawisko w tym przypadku?
5. Wartość ROC przypisz do zmiennej auroc.

Uwaga:

• koniecznie ustaw verbose=-1 przy tworzeniu LGBMClassifier, żeby uniknąć kolosalnej ilości logów, która potrafi też wyłączyć Jupytera
• pamiętaj o ustawieniu importance_type, random_state=0 i n_jobs, oraz ewentualnie verbose w RandomizedSearchCV dla śledzenia przebiegu
• istnieje możliwość, że ustawienie n_jobs dla grid searcha będzie szybsze niż dla samego LightGBM; odpowiada to tuningowi wielu klasyfikatorów równolegle, przy wolniejszym treningu pojedynczych klasyfikatorów
• nie ustawiaj wszędzie n_jobs=-1, bo wtedy stworzysz więcej procesów niż rdzeni i spowodujesz thread contention

from sklearn.metrics import classification_report


# function for printing the classification report
def get_classification_report(model, X_test, y_test, title):

    y_pred = model.predict(X_test)

    print(title)

    print(classification_report(y_true=y_test, y_pred=y_pred))

get_classification_report(
    LGBM_clf,
    X_test,
    y_test,
    "Classification report for LGBM without hyper-parameter tuning:",
)

Classification report for LGBM without hyper-parameter tuning:
              precision    recall  f1-score   support

           0       0.98      0.98      0.98      2002
           1       0.60      0.60      0.60        99

    accuracy                           0.96      2101
   macro avg       0.79      0.79      0.79      2101
weighted avg       0.96      0.96      0.96      2101

# your_code

from sklearn.model_selection import RandomizedSearchCV

param_grid = {
    "n_estimators": [100, 250, 500],
    "learning_rate": [0.05, 0.1, 0.2],
    "num_leaves": [31, 48, 64],
    "colsample_bytree": [0.8, 0.9, 1.0],
    "subsample": [0.8, 0.9, 1.0],
}

LGBM_estimator = LGBMClassifier(
    importance_type="gain", random_state=0, n_jobs=-1, verbose=-1
)

LGBM_RS = RandomizedSearchCV(
    estimator=LGBM_estimator,
    cv=5,
    param_distributions=param_grid,
    verbose=0,
    scoring="roc_auc",
    n_iter=30,
)

LGBM_RS.fit(X_reb, y_reb)

RandomizedSearchCV(cv=5,
                   estimator=LGBMClassifier(importance_type='gain', n_jobs=-1,
                                            random_state=0, verbose=-1),
                   n_iter=30,
                   param_distributions={'colsample_bytree': [0.8, 0.9, 1.0],
                                        'learning_rate': [0.05, 0.1, 0.2],
                                        'n_estimators': [100, 250, 500],
                                        'num_leaves': [31, 48, 64],
                                        'subsample': [0.8, 0.9, 1.0]},
                   scoring='roc_auc')

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

print(f"Best parameters: {LGBM_RS.best_params_}\n")

get_classification_report(
    LGBM_RS,
    X_test,
    y_test,
    "Classification report for LGBM with hyper-parameter tuning:",
)

Best parameters: {'subsample': 1.0, 'num_leaves': 48, 'n_estimators': 500, 'learning_rate': 0.2, 'colsample_bytree': 0.9}

Classification report for LGBM with hyper-parameter tuning:
              precision    recall  f1-score   support

           0       0.98      0.99      0.99      2002
           1       0.80      0.54      0.64        99

    accuracy                           0.97      2101
   macro avg       0.89      0.76      0.81      2101
weighted avg       0.97      0.97      0.97      2101

auroc = get_auroc_score(LGBM_RS, X_test, y_test)

assert 0.9 <= auroc <= 0.99

AUROC score: 0.9397925307016216

// skomentuj tutaj

To czy jest to pożądane zjawisko, to zależy od naszej taktyki. Załóżmy, że bardziej nas interesuje precyzja i chcemy być maksylmanie skuteczni w naszych przewidywaniach bankructwa, bo fałszywy alarm może prowadzić np. do zamknięcia dobrze prosperującej firmy albo drastycznej restrukturyzacji. W takim wypadku, dostrojenie hiper-parametrów miało duże znaczenie. Precyzja znacznie się polepszyła, co jest dla nas bardzo korzystne. Warto zwrócić uwagę na to, że AUROC w obu przypadkach jest dosyć podobny (~0.94). Podczas wykonywania tego labolatorium napotkałem masę różnych artykułów pt. "Why you should not use AUROC". Autorzy słusznie zwracali uwagę na potrzebę oceny modelu również innymi metrykami (często przewijał się Precision-Recall Curve), co w tym wypadku ma swoje odzwierciedlenie, jako że dla dwóch bliskich wyników AUROC, inne metryki, z perspektywy naszego biznesplanu - ważniejsze, mocno się różniły.

Boosting - podsumowanie

1. Model oparty o uczenie zespołowe
2. Kolejne modele są dodawane sekwencyjnie i uczą się na błędach poprzedników
3. Nauka typowo jest oparta o minimalizację funkcji kosztu (błędu), z użyciem spadku wzdłuż gradientu
4. Wiodący model klasyfikacji dla danych tabelarycznych, z 2 głównymi implementacjami: XGBoost i LightGBM
5. Liczne hiperparametry, wymagające odpowiednich metod dostrajania

Wyjaśnialna AI

W ostatnich latach zaczęto zwracać coraz większą uwagę na wpływ sztucznej inteligencji na społeczeństwo, a na niektórych czołowych konferencjach ML nawet obowiązkowa jest sekcja "Social impact" w artykułach naukowych. Typowo im lepszy model, tym bardziej złożony, a najpopularniejsze modele boostingu są z natury skomplikowane. Kiedy mają podejmować krytyczne decyzje, to musimy wiedzieć, czemu predykcja jest taka, a nie inna. Jest to poddziedzina uczenia maszynowego - wyjaśnialna AI (explainable AI, XAI).

Taka informacja jest cenna, bo dzięki temu lepiej wiemy, co robi model. Jest to ważne z kilku powodów:

1. Wymogi prawne - wdrażanie algorytmów w ekonomii, prawie etc. ma coraz częściej konkretne wymagania prawne co do wyjaśnialności predykcji
2. Dodatkowa wiedza dla użytkowników - często dodatkowe obserwacje co do próbek są ciekawe same w sobie i dają wiedzę użytkownikowi (często posiadającemu specjalistyczną wiedzę z dziedziny), czasem nawet bardziej niż sam model predykcyjny
3. Analiza modelu - dodatkowa wiedza o wewnętrznym działaniu algorytmu pozwala go lepiej zrozumieć i ulepszyć wyniki, np. przez lepszy preprocessing danych

W szczególności można ją podzielić na globalną oraz lokalną interpretowalność (global / local interpretability). Ta pierwsza próbuje wyjaśnić, czemu ogólnie model działa tak, jak działa. Analizuje strukturę modelu oraz trendy w jego predykcjach, aby podsumować w prostszy sposób jego tok myślenia. Interpretowalność lokalna z kolei dotyczy predykcji dla konkretnych próbek - czemu dla danego przykładu model podejmuje dla niego taką, a nie inną decyzję o klasyfikacji.

W szczególności podstawowym sposobem interpretowalności jest ważność cech (feature importance). Wyznacza ona, jak ważne są poszczególne cechy:

• w wariancie globalnym, jak mocno model opiera się na poszczególnych cechach
• w wariancie lokalnym, jak mocno konkretne wartości cech wpłynęły na predykcję, i w jaki sposób

Teraz będzie nas interesować globalna ważność cech. Dla modeli drzewiastych definiuje się ją bardzo prosto. Każdy podział w drzewie decyzyjnym wykorzystuje jakąś cechę, i redukuje z pomocą podziału funkcję kosztu (np. entropię) o określoną ilość. Dla drzewa decyzyjnego ważność to sumaryczna redukcja entropii, jaką udało się uzyskać za pomocą danej cechy. Dla lasów losowych i boostingu sumujemy te wartości dla wszystkich drzew. Alternatywnie można też użyć liczby splitów, w jakiej została użyta dana cecha, ale jest to mniej standardowe.

Warto zauważyć, że taka ważność cech jest względna:

• nie mówimy, jak bardzo ogólnie ważna jest jakaś cecha, tylko jak bardzo przydatna była dla naszego modelu w celu jego wytrenowania
• ważność cech można tylko porównywać ze sobą, np. jedna jest 2 razy ważniejsza od drugiej; nie ma ogólnych progów ważności

Ze względu na powyższe, ważności cech normalizuje się często do zakresu [0, 1] dla łatwiejszego porównywania.

Zadanie 9 (0.5 punktu)

1. Wybierz 5 najważniejszych cech dla drzewa decyzyjnego. Przedstaw wyniki na poziomym wykresie słupkowym. Użyj czytelnych nazw cech ze zmiennej feature_names.
2. Powtórz powyższe dla lasu losowego, oraz dla boostingu (tutaj znormalizuj wyniki - patrz uwaga niżej). Wybierz te hiperparametry, które dały wcześniej najlepsze wyniki.
3. Skomentuj, czy wybrane cechy twoim zdaniem mają sens jako najważniejsze cechy.

Uwaga: Scikit-learn normalizuje ważności do zakresu [0, 1], natomiast LightGBM nie. Musisz to znormalizować samodzielnie, dzieląc przez sumę.

# function for plotting top 5 feature importances on a specific ax
def plot_feature_importances(feature_importances, ax, title):

    head_indices = np.argsort(feature_importances)[-5:]
    head = feature_importances[head_indices]
    head_names = feature_names[head_indices]

    ax.bar(head_names, head)
    ax.set_ylabel("Feature Importance")
    ax.set_title(title)
    ax.tick_params(rotation=70, axis="x")

# your_code
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(30, 9))

tree_feature_importances = decision_tree.feature_importances_

forest_feature_importances = random_forest.feature_importances_

LGBM_feature_importances = LGBM_clf.feature_importances_
LGBM_feature_importances /= LGBM_feature_importances.sum()

plot_feature_importances(
    tree_feature_importances, ax1, "Feature Importance from Decision Tree [Top 5]"
)
plot_feature_importances(
    forest_feature_importances, ax2, "Feature Importance from Random Forest [Top 5]"
)
plot_feature_importances(
    LGBM_feature_importances, ax3, "Feature Importance from LGBM [Top 5]"
)

// skomentuj tutaj

W trzech różnych metodach zostały wyróżnione całkiem "logiczne" cechy: takie jak sales(n)/sales(n-1) oraz profit on operating activites/financial expenses. To ważne cechy, których odpowiednie stosunki (tzn. $\geq 1$) świadczą o dobrej kondycji firmy. Pierwsza kładzie nacisk na zdolność do rozwoju, a druga na wydajność finansową firmy (potrzebną np. w przypadku długów). Myslę, że wybór tych cech jak najbardziej ma sens.

Dla zainteresowanych

Najpopularniejszym podejściem do interpretowalności lokalnych jest SHAP (SHapley Additive exPlanations), metoda oparta o kooperatywną teorię gier. Traktuje się cechy modelu jak zbiór graczy, podzielonych na dwie drużyny (koalicje): jedna chce zaklasyfikować próbkę jako negatywną, a druga jako pozytywną. O ostatecznej decyzji decyduje model, który wykorzystuje te wartości cech. Powstaje pytanie - w jakim stopniu wartości cech przyczyniły się do wyniku swojej drużyny? Można to obliczyć jako wartości Shapleya (Shapley values), które dla modeli ML oblicza algorytm SHAP. Ma on bardzo znaczące, udowodnione matematycznie zalety, a dodatkowo posiada wyjątkowo efektywną implementację dla modeli drzewiastych oraz dobre wizualizacje.

Bardzo intuicyjnie, na prostym przykładzie, SHAPa wyjaśnia pierwsza część tego artykułu. Dobrze i dość szczegółówo SHAPa wyjaśnia jego autor w tym filmie.

Wyjaśnialna AI - podsumowanie

1. Problem zrozumienia, jak wnioskuje model i czemu podejmuje dane decyzje
2. Ważne zarówno z perspektywy data scientist'a, jak i użytkowników systemu
3. Można wyjaśniać model lokalnie (konkretne predykcje) lub globalnie (wpływ poszczególnych cech)

Zadanie 10 dla chętnych (3 punkty)

Dokonaj selekcji cech, usuwając 20% najsłabszych cech. Może się tu przydać klasa SelectPercentile. Czy Random Forest i LightGBM (bez dostrajania hiperparametrów, dla uproszczenia) wytrenowane bez najsłabszych cech dają lepszy wynik (AUROC lub innej metryki)?

Wykorzystaj po 1 algorytmie z 3 grup algorytmów selekcji cech:

1. Filter methods - mierzymy ważność każdej cechy niezależnie, za pomocą pewnej miary (typowo ze statystyki lub teorii informacji), a potem odrzucamy (filtrujemy) te o najniższej ważności. Są to np. chi2 i mutual_info_classif z pakietu sklearn.feature_selection.
2. Embedded methods - klasyfikator sam zwraca ważność cech, jest jego wbudowaną cechą (stąd nazwa). Jest to w szczególności właściwość wszystkich zespołowych klasyfikatorów drzewiastych. Mają po wytrenowaniu atrybut feature_importances_.
3. Wrapper methods - algorytmy wykorzystujące w środku używany model (stąd nazwa), mierzące ważność cech za pomocą ich wpływu na jakość klasyfikatora. Jest to np. recursive feature elimination (klasa RFE). W tym algorytmie trenujemy klasyfikator na wszystkich cechach, wyrzucamy najsłabszą, trenujemy znowu i tak dalej.

Typowo metody filter są najszybsze, ale dają najsłabszy wynik, natomiast metody wrapper są najwolniejsze i dają najlepszy wynik. Metody embedded są gdzieś pośrodku.

Dla zainteresowanych, inne znane i bardzo dobre algorytmy:

• Relief (filter method) oraz warianty, szczególnie ReliefF, SURF i MultiSURF (biblioteka ReBATE): Wikipedia), artykuł "Benchmarking Relief-Based Feature Selection Methods"
• Boruta (wrapper method), stworzony na Uniwersytecie Warszawskim, łączący Random Forest oraz testy statystyczne (biblioteka boruta_py): link 1, link 2

from sklearn.feature_selection import RFE, SelectPercentile, mutual_info_classif

print(f"Without data selection\n")

print(f"Random forest:")

random_forest = RandomForestClassifier(
    n_estimators=500, criterion="entropy", random_state=0, n_jobs=-1
)

random_forest.fit(X_reb, y_reb)

get_auroc_score(random_forest, X_test, y_test)
get_classification_report(
    random_forest, X_test, y_test, "Classification report for Random forest:"
)

print(f"LGBM:")

LGBM = LGBMClassifier(importance_type="gain", random_state=0, n_jobs=-1)

LGBM.fit(X_reb, y_reb)

get_auroc_score(LGBM, X_test, y_test)
get_classification_report(LGBM, X_test, y_test, "Classification report for LGBM:")

Without data selection

Random forest:
AUROC score: 0.9047644274917003
Classification report for Random forest:
              precision    recall  f1-score   support

           0       0.97      0.98      0.98      2002
           1       0.53      0.40      0.46        99

    accuracy                           0.96      2101
   macro avg       0.75      0.69      0.72      2101
weighted avg       0.95      0.96      0.95      2101

LGBM:
AUROC score: 0.9433748070111706
Classification report for LGBM:
              precision    recall  f1-score   support

           0       0.98      0.98      0.98      2002
           1       0.60      0.60      0.60        99

    accuracy                           0.96      2101
   macro avg       0.79      0.79      0.79      2101
weighted avg       0.96      0.96      0.96      2101

print(f"Filter Methods\n")

filter_method = SelectPercentile(score_func=mutual_info_classif, percentile=80)

X_train_selected = filter_method.fit_transform(X_reb, y_reb)
X_test_selected = filter_method.transform(X_test)

print(f"Random forest:")

random_forest = RandomForestClassifier(
    n_estimators=500, criterion="entropy", random_state=0, n_jobs=-1
)

random_forest.fit(X_train_selected, y_reb)

get_auroc_score(random_forest, X_test_selected, y_test)
get_classification_report(
    random_forest, X_test_selected, y_test, "Classification report for Random forest:"
)

print(f"LGBM:")

LGBM = LGBMClassifier(importance_type="gain", random_state=0, n_jobs=-1)

LGBM.fit(X_train_selected, y_reb)

get_auroc_score(LGBM, X_test_selected, y_test)
get_classification_report(
    LGBM, X_test_selected, y_test, "Classification report for LGBM:"
)

Filter Methods

Random forest:
AUROC score: 0.8951906679179406
Classification report for Random forest:
              precision    recall  f1-score   support

           0       0.97      0.98      0.97      2002
           1       0.46      0.42      0.44        99

    accuracy                           0.95      2101
   macro avg       0.71      0.70      0.71      2101
weighted avg       0.95      0.95      0.95      2101

LGBM:
AUROC score: 0.9245602881966518
Classification report for LGBM:
              precision    recall  f1-score   support

           0       0.98      0.97      0.98      2002
           1       0.50      0.55      0.52        99

    accuracy                           0.95      2101
   macro avg       0.74      0.76      0.75      2101
weighted avg       0.95      0.95      0.95      2101

# function for selecting 80% of the most importatnt features
def get_selected_features(feature_importances, X_train, X_test):
    # get top 80% feature count
    count = int(feature_importances.size * 0.8)

    # calculate the threshold and cut the least important 20%
    threshold = np.min(np.sort(feature_importances)[count:])
    selected = feature_importances >= threshold
    X_train_selected = X_train[:, selected]
    X_test_selected = X_test[:, selected]

    return X_train_selected, X_test_selected

print(f"Embedded Methods\n")

print(f"Random forest:")

X_train_forest, X_test_forest = get_selected_features(
    forest_feature_importances, X_reb, X_test
)

random_forest = RandomForestClassifier(
    n_estimators=500, criterion="entropy", random_state=0, n_jobs=-1
)

random_forest.fit(X_train_forest, y_reb)

get_auroc_score(random_forest, X_test_forest, y_test)
get_classification_report(
    random_forest, X_test_forest, y_test, "Classification report for Random forest:"
)

print(f"LGBM:")

X_train_LGBM, X_test_LGBM = get_selected_features(
    LGBM_feature_importances, X_reb, X_test
)

LGBM = LGBMClassifier(importance_type="gain", random_state=0, n_jobs=-1)

LGBM.fit(X_train_LGBM, y_reb)

get_auroc_score(LGBM, X_test_LGBM, y_test)
get_classification_report(LGBM, X_test_LGBM, y_test, "Classification report for LGBM:")

Embedded Methods

Random forest:
AUROC score: 0.9151757333575514
Classification report for Random forest:
              precision    recall  f1-score   support

           0       0.98      0.97      0.97      2002
           1       0.44      0.54      0.48        99

    accuracy                           0.95      2101
   macro avg       0.71      0.75      0.73      2101
weighted avg       0.95      0.95      0.95      2101

LGBM:
AUROC score: 0.9339397975761612
Classification report for LGBM:
              precision    recall  f1-score   support

           0       0.98      0.96      0.97      2002
           1       0.45      0.62      0.52        99

    accuracy                           0.95      2101
   macro avg       0.71      0.79      0.75      2101
weighted avg       0.96      0.95      0.95      2101

print(f"Wrapper Methods\n")

print(f"Random forest:")
random_forest = RandomForestClassifier(
    n_estimators=500, criterion="entropy", random_state=0, n_jobs=-1
)

random_forest_RFE = RFE(estimator=random_forest, n_features_to_select=0.8)

random_forest_RFE.fit(X_reb, y_reb)

get_auroc_score(random_forest_RFE, X_test, y_test)
get_classification_report(
    random_forest_RFE, X_test, y_test, "Classification report for Random forest:"
)

print(f"LGBM:")
LGBM = LGBMClassifier(importance_type="gain", random_state=0, n_jobs=-1)

LGBM_RFE = RFE(estimator=LGBM, n_features_to_select=0.8)

LGBM_RFE.fit(X_reb, y_reb)

get_auroc_score(LGBM_RFE, X_test, y_test)
get_classification_report(LGBM_RFE, X_test, y_test, "Classification report for LGBM:")

Wrapper Methods

Random forest:
AUROC score: 0.9105187741551377
Classification report for Random forest:
              precision    recall  f1-score   support

           0       0.97      0.98      0.98      2002
           1       0.51      0.41      0.46        99

    accuracy                           0.95      2101
   macro avg       0.74      0.70      0.72      2101
weighted avg       0.95      0.95      0.95      2101

LGBM:
AUROC score: 0.9436876255058073
Classification report for LGBM:
              precision    recall  f1-score   support

           0       0.98      0.98      0.98      2002
           1       0.61      0.58      0.59        99

    accuracy                           0.96      2101
   macro avg       0.79      0.78      0.79      2101
weighted avg       0.96      0.96      0.96      2101

Zgodnie z założeniami metody filter okazały się najsłabsze i tylko pogorszyły model. Metody Wrapper najbardziej usprawniły model, a metody embedded były gdzieś pośrodku. W przypadku innych metryk: precyzja i pokrycie, selekcja cech na ogół negatywnie na nie wpłynęła, bądź lekko polepszyła.